Carotenoid deposition simulation for common yellowthroat feather coloration.
Conor C. Taff
Purpose
The goal here is to demonstrate how different patterns of variation in i) the shape of structural white spectra, ii) the UV shift of the white spectra, iii) amount of carotenoid deposition, and iv) the overall reflectance level generate different patterns of covariance in downstream metrics of coloration. A simulation is set up allowing me to vary each of these characteristics separately or together to demonstrate the effects.
Simulation
The simulation here is simple and involves four stages that represent the four possible sources of variation in ultimate color measurements. These four transformations are employed sequentially in the order presented to a set number of randomly drawn individuals. The amount of variation in each process can be set to examine the consequences of variation in each process alone or any combination of processes.
1. Shape of spectral curve
To generate the structural shape of the white feather spectra, I start with the average white reflectance curve in Shawkey & Hill (Shawkey & Hill, 2005, fig. 1A ‘before treatment’ and ‘after treatment’ can be chosen). I extracted these curves using WebPlotDigitizer and interpolated points for every 1 nm wavelength by fitting a smoothed regression. This results in the ‘typical’ white color that is used as the default when structural color does not vary (figure @ref(fig:white-shawkey)).
# load shawkey figure data
shaw <- read.delim(here::here("1_raw_data", "shawkey_2005_fig1.txt"))
shaw_a <- subset(shaw, shaw$type == "after")
shaw_b <- subset(shaw, shaw$type == "before")
# Fit loess smooth and interpolate value at each 1 nm
r <- data.frame(wavelength = seq(300, 700, 1))
sa_lo <- loess(reflectance ~ wavelength, data = shaw_a, span = 0.1, surface = "direct")
sa_pred <- data.frame(wavelength = r$wavelength, reflectance = predict(sa_lo, newdata = r))
sb_lo <- loess(reflectance ~ wavelength, data = shaw_b, span = 0.2, surface = "direct")
sb_pred <- data.frame(wavelength = r$wavelength, reflectance = predict(sb_lo, newdata = r))
s_pred <- data.frame(group = rep(c("Before", "After"), each = nrow(r)),
wavelength = rep(r$wavelength, 2),
reflectance = c(sb_pred$reflectance, sa_pred$reflectance))
# Plot the lines
ggplot(data = s_pred, mapping = aes(x = wavelength, y = reflectance, color = group)) +
geom_smooth(se = FALSE, span = 0.1) +
scale_color_manual(values = c("orange", "slateblue")) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = c(0.8, 0.2)) +
theme(legend.title = element_blank()) +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)")
Structural white curves taken from Shawkey and Hill 2005 Figure 1A before and after feather bleaching treatments.
To generate individuals that differ in their underlying structural shape, I multiplied this entire spectra by a percentage chosen from a random normal distribution with mean of 100% and SD of 15%. Multiplying by a percentage has the effect of both raising or lowering the overall curve and also of changing the overall percentage of reflectance that is represented in different wavelengths. The different spectral shapes generated by this sampling process can bee seen in Figure @ref(fig:white-variation).
# multiply by percentages
pop_p <- rnorm(25, 1, 0.15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
pop_r[, i] <- sb_pred$reflectance * rep(pop_p[i], nrow(sb_pred))
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p1 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.1) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("Before") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 8)
# multiply by percentages
pop_p <- rnorm(25, 1, 0.15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
pop_r[, i] <- sa_pred$reflectance * rep(pop_p[i], nrow(sb_pred))
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p2 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.1) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("After") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 8)
ggpubr::ggarrange(p2, p1)
Sampled spectral shapes starting from the base curves from Shawkey and Hill 2005 and multiplying by a random normal sample of percentages with mean of 100% and SD of 15%. Panel A starts with the base curve for ‘after bleaching’ and panel B starts with the base curve for ‘before bleaching.’
Amount of UV in spectral curve
I’m not sure this is really a distinct mechanism from the shape section above, but we also discussed the fact that the same curve could be shifted to have more UV reflectance. I’ve implemented this source of variation in the following way. First, for each wavelength across the spectra I calculate one minus the percentage of the maximum reflectance for that wavelength (so low reflectance areas in UV end up with high numbers). Second, I draw a number from a random normal distribution with mean of 0 (no distortion) and SD of 0.15. I multiply this fixed value by the vector calculated above and add it to the original spectra. This results in a shifted towards UV peak for positive draws and a shifted away from UV peak for negative draws (Figure @ref(fig:uv-variation)).
For now, in the the simulations below I am ignoring this mechanism for simplicity.
# multiply by uv shift
pop_p <- rnorm(25, 0, 0.15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
temp <- (1 - sb_pred$reflectance / rep(max(sb_pred$reflectance), nrow(sb_pred)))
draw <- rnorm(1, 0, 12)
new <- sb_pred$reflectance + rep(draw, nrow(sb_pred)) * temp
new <- pmax(new, 0)
pop_r[, i] <- new
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p1 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.05) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("Before") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 8)
# multiply by uv shift
pop_p <- rnorm(25, 0, 0.15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
temp <- (1 - sa_pred$reflectance / rep(max(sa_pred$reflectance), nrow(sa_pred)))
draw <- rnorm(1, 0, 12)
new <- sa_pred$reflectance + rep(draw, nrow(sa_pred)) * temp
new <- pmax(new, 0)
pop_r[, i] <- new
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p2 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.05) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("After") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 8)
ggpubr::ggarrange(p2, p1) 
Sampled spectral shapes with shifts indicating more or less relative UV while max reflectance stays the same. Panels show the process applied to the ‘after’ (A) or ‘before’ (B) spectra from Shawkey and Hill 2005.
3. Overall level of reflectance
Here the goal was to change the overall level of reflectance while maintaining the shape of the underlying spectral shape. To accomplish this, I added or subtracted a fixed percentage of reflectance to the entire curve. When generating between-individual differences in this component, I chose a fixed value to add or subtract from a random normal distribution with mean of 0 and SD of 15%. In some cases, this subtraction could lead to negative reflectance values, so after the transformation I set any reflectance values less than zero to zero. This can result in slightly different shapes to curves at the low end, but is necessary to create realistic curves. Clipping can be avoided by starting with a curve that has a higher initial reflectance or changing the random addition numbers (e.g., use a mean of 10 and SD of 10%). Holding everything else constant, this process results in a population of curves that look like Figure @ref(fig:reflect-variation).
# add random amount
pop_p <- rnorm(25, 0, 15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
pop_r[, i] <- sb_pred$reflectance + rep(pop_p[i], nrow(sb_pred))
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p1 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.1) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("Before") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 8)
# add percentages
pop_p <- rnorm(25, 0, 15)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(s_pred)))
for(i in 1:length(pop_p)){
pop_r[, i] <- sa_pred$reflectance + rep(pop_p[i], nrow(sb_pred))
}
pop_r$wavelength <- s_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p2 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
geom_smooth(se = FALSE, span = 0.1) +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("After") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 8)
ggpubr::ggarrange(p2, p1) 
Sampled spectral shapes starting from the base curves from Shawkey and Hill 2005 and adding a random normal sample of percentages with mean of 0 and SD of 15%. Panel A starts with the base curve for ‘after bleaching’ and panel B starts with the base curve for ‘before bleaching.’
4. Amount of carotenoid deposition
Last, we can take the curves generated by the two processes above and then layer carotenoid deposition on top of them. In common yellowthroats, yellow color is generated solely by lutein. I found a plot of the absorption spectrum of lutein online (there are many published spectra) and digitized it using WebPlotDigitizer as described above so that I ended up with an absorption value for every wavelength from 300 to 700 nm (Figure @ref(fig:lutein-only)).
# Read in lutein data
# load lutein data that I copied from an absorbance spectra posted online
lut <- read.delim(here::here("1_raw_data", "lutein.txt"))
# In order to make it flat outside of range I'm adding more points that extend out at the beginning and end
lut_begin <- data.frame(molecule = "lutein",
wavelength = seq(280, 310, 5),
absorbance = 0.05209852)
lut_end <- data.frame(molecule = "lutein",
wavelength = seq(565, 800, 5),
absorbance = 0.008027346)
lut <- rbind(lut, lut_begin, lut_end)
# Fit loess smooth and interpolate value at each 1 nm
r <- data.frame(wavelength = seq(300, 700, 1))
lu_lo <- loess(absorbance ~ wavelength, data = lut, span = 0.05, surface = "direct")
lu_pred <- data.frame(wavelength = r$wavelength, absorbance = predict(lu_lo, newdata = r))
# Plot the lines
ggplot(data = lu_pred, mapping = aes(x = wavelength, y = absorbance)) +
#geom_smooth(se = FALSE, span = 0.04, color = "goldenrod", size = 1.5) +
#geom_point(alpha = 0.5) +
geom_line(color = "goldenrod", size = 1.5) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = c(0.8, 0.2)) +
theme(legend.title = element_blank()) +
xlab("Wavelength (nm)") +
ylab("Absorbance (%)")
Absorbance spectrum of lutein.
The absorbance spectrum that I used was calibrated such that the maximum absorbance at 452 nm was ~96%. I generated samples of different amounts of lutein deposition by multiplying this lutein spectrum by number chosen from a random normal distribution with mean of 0.8 and SD of 0.1. Spectral curves with carotenoid absorption were then calculated by the formula: white - white * lutein absorption. When underlying white coloration and reflectance is held constant, this process resulted in a population of spectral shapes as shown in Figure @ref(fig:lutein-variation). Note that there is no assumption here about the actual amount of carotenoid deposition represented by these curves; rather, it is designed to create realistic spectral shapes and variation without considering the detailed physical properties. As above, I clipped any negative reflectance values to 0.
# multiply by lutein
pop_p <- rnorm(25, 0.8, 0.1)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(sb_pred)))
for(i in 1:length(pop_p)){
draw <- rnorm(1, 0.7, 0.15)
if(draw > 1){draw <- 1}
temp <- draw * lu_pred$absorbance
pop_r[, i] <- sb_pred$reflectance - (sb_pred$reflectance * temp)
}
pop_r$wavelength <- sb_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p1 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
#geom_smooth(se = FALSE, span = 0.05) +
geom_line() +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("Before") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "B", size = 8)
# multiply by lutein
pop_p <- rnorm(25, 0.8, 0.1)
pop_r <- as.data.frame(matrix(ncol = length(pop_p), nrow = nrow(sa_pred)))
for(i in 1:length(pop_p)){
draw <- rnorm(1, 0.7, 0.15)
if(draw > 1){draw <- 1}
temp <- draw * lu_pred$absorbance
pop_r[, i] <- sa_pred$reflectance - (sa_pred$reflectance * temp)
}
pop_r$wavelength <- sa_pred$wavelength
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
# Plot the lines
p2 <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
#geom_smooth(se = FALSE, span = 0.05) +
geom_line() +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)") +
ggtitle("After") +
annotate(geom = "text", x = -Inf, y = Inf, hjust = -.5, vjust = 1.5, label = "A", size = 8)
ggpubr::ggarrange(p2, p1) 
Variation in lutein deposition added to the basic structure for ‘after’ (A) or before (B) white spectra.
5. Measurement error
This just raises or lowers the entire curve shape by a random amount after all the prior transformations are finished.
Calculating color metrics
The three step process described above can generate an arbitrary number of spectra that vary in any way desired (i.e., we can specify variation in any one, two, or all three of the contributing processes). After generating a sample of spectra as desired, I calculated a variety of downstream color metrics that summarize certain properties of the color spectra. For now I’ve calculated metrics exactly as we describe in Freeman-Gallant et al. ((2010)). This includes total brightness, uv brightness, yellow brightness, uv saturation, yellow saturation, and Ccar. The final result is a table of coloration metrics for a population of a defined sample size and set sources of variation in overall color that can be used to ask how each metric covaries.
Function
The description above describes how each step in the simulation works, but the entire set of steps is contained in a single function called lutein_sim. The function takes arguments for the sample size desired, which initial curve to start with, and the mean and standard deviation settings for each source of variation described above. The output includes i) the full spectra for each individual, ii) the summary table of color metrics, iii) a plot with the full spectra for each individual, and iv) a correlation plot of the key pair-wise color metrics. Each of these four outputs can be accessed from the list object returned by the function. The entire set of steps describe above is run in a single line of code. Because this is drawing from a random distribution, remember to set a seed if you want identical results to be produced each time. For now I just have the function running as a standalone script, but I can set it up to install via GitHub if we want. You’ll need to have the following packages already installed in your R environment: tidyverse; ggplot2; viridis; GGally. Once those are installed, you can just run the single code block from the repository to load the function and then have access to the function to generate new variations.
Note that I bumped up the overall reflectance a bit from the Shawkey figure just to stop it from clipping at 0 when adding variation. That should not change any correlation patterns.
ADDITION 7/26. I’ve cleaned up the code a bit and added the ability to specify correlations between each of the mechanisms described above. By default correlations are set to 0 (same as previous version), but each can now be specified individually. Some of the argument names are changed for clarity.
carotenoid_sim <- function(n = 50, shape1_mu = 1, shape1_sd = 0,
shape2_mu = 0, shape2_sd = 0,
totalb_mu = 12, totalb_sd = 0,
pigment_mu = 0.75, pigment_sd = 0,
shape1_shape2_cor = 0, shape1_totalb_cor = 0, shape1_pigment_cor = 0,
shape2_totalb_cor = 0, shape2_pigment_cor = 0, totalb_pigment_cor = 0,
base_shape = "after", measure_error = 0){
# Convert uv shift (shape 2) and total brightness sds to *100 since the percentages are on whole number scales
shape2_sd <- shape2_sd * 100
totalb_sd <- totalb_sd * 100
# Build a variance covariance matrix based on the numbers listed above
sim_matrix <- matrix(nrow = 4, ncol = 4, data = c(
shape1_sd^2, shape1_shape2_cor*shape1_sd*shape2_sd, shape1_totalb_cor*shape1_sd*totalb_sd, shape1_pigment_cor*shape1_sd*pigment_sd,
shape1_shape2_cor*shape1_sd*shape2_sd, shape2_sd^2, shape2_totalb_cor*shape2_sd*totalb_sd, shape2_pigment_cor*shape2_sd*pigment_sd,
shape1_totalb_cor*shape1_sd*totalb_sd, shape2_totalb_cor*shape2_sd*totalb_sd, totalb_sd^2, totalb_pigment_cor*totalb_sd*pigment_sd,
shape1_pigment_cor*shape1_sd*pigment_sd, shape2_pigment_cor*shape2_sd*pigment_sd, totalb_pigment_cor*totalb_sd*pigment_sd, pigment_sd^2
))
# Sample values for individuals from a multivariate normal distribution based on the var-cov matrix constructed above
individuals <- as.data.frame(mvrnorm(n = n, mu = c(shape1_mu, shape2_mu, totalb_mu, pigment_mu), Sigma = sim_matrix))
colnames(individuals) <- c("shape1", "shape2", "totalb", "pigment")
# Pigment values above 1 result in negative reflectance. Clip these out. Changing mean/sd settings in call can prevent these from happening
individuals$pigment <- ifelse(individuals$pigment > 1, 1, individuals$pigment)
# Set up for using before or after spectra
if(base_shape == "after"){bs <- "after"}
if(base_shape == "before"){bs <- "before"}
# Read in basic shape data. This requires the file to be present. Should really be hardcoded in.
# load shawkey figure data
shaw <- read.delim(here::here("1_raw_data", "shawkey_2005_fig1.txt"))
shaw <- subset(shaw, shaw$type == bs)
# Fit and predict from loess smoothed line to get reflectance at each nm wavelength
r <- data.frame(wavelength = seq(300, 700, 1))
lo_fit <- loess(reflectance ~ wavelength, data = shaw, span = 0.1, surface = "direct")
lo_pred <- data.frame(wavelength = r$wavelength, reflectance = predict(lo_fit, newdata = r))
# Read in lutein shape data. This also requires file and should be hardcoded in. I should add different pigments.
# load lutein data that I copied from an absorbance spectra posted online
lut <- read.delim(here::here("1_raw_data", "lutein.txt"))
# In order to make it flat outside of range I'm adding more points that extend out at the beginning and end
lut_begin <- data.frame(molecule = "lutein",
wavelength = seq(280, 310, 5),
absorbance = 0.05209852)
lut_end <- data.frame(molecule = "lutein",
wavelength = seq(565, 800, 5),
absorbance = 0.008027346)
lut <- rbind(lut, lut_begin, lut_end)
# Fit loess smooth and interpolate value at each 1 nm
lu_lo <- loess(absorbance ~ wavelength, data = lut, span = 0.05, surface = "direct")
lu_pred <- data.frame(wavelength = r$wavelength, absorbance = predict(lu_lo, newdata = r))
# 1. Add in white variation
# multiply by percentages
pop_r <- as.data.frame(matrix(ncol = nrow(individuals), nrow = nrow(lo_pred)))
for(i in 1:n){
pop_r[, i] <- lo_pred$reflectance * rep(individuals[i, 1], nrow(lo_pred))
}
# 2. Add in uv shift variation
# multiply by uv shift
for(i in 1:n){
temp <- (1 - pop_r[, i] / rep(max(pop_r[, i]), nrow(lo_pred)))
draw <- individuals[i, 2]
new <- pop_r[, i] + rep(draw, nrow(lo_pred)) * temp
pop_r[, i] <- new
}
# 3. Add in overall reflectance variation
pop_p2 <- individuals$totalb
for(i in 1:n){
pop_r[, i] <- pop_r[, i] + rep(pop_p2[i], nrow(lo_pred))
}
# 4. Add in lutein deposition variation
# multiply by lutein
for(i in 1:n){
draw <- individuals$pigment[i]
temp <- draw * lu_pred$absorbance
pop_r[, i] <- pop_r[, i] - (pop_r[, i] * temp)
}
# Add measurement error
measure_error <- measure_error*100 # convert to percentage scale
errors <- rnorm(n, mean = 0, sd = measure_error)
for(i in 1:n){
pop_r[, i] <- pop_r[, i] + rep(errors[i], nrow(pop_r))
}
pop_r[pop_r < 0] <- 0
# Add back in a column with wavelength information
pop_r$wavelength <- lo_pred$wavelength
pop_r[pop_r < 0] <- 0
pop_r2 <- pivot_longer(pop_r, cols = starts_with("V"), names_to = "id", values_to = "reflectance")
## Make figure
p_spec <- ggplot(data = pop_r2, mapping = aes(x = wavelength, y = reflectance, color = id)) +
#geom_smooth(se = FALSE, span = 0.05) +
geom_line() +
scale_color_viridis(discrete = TRUE) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
legend.title = element_text(size = 14)) +
guides(color = "none") +
xlab("Wavelength (nm)") +
ylab("Reflectance (%)")
## make dataframe of color metrics
color_met <- data.frame(id = unique(pop_r2$id),
b_tot = rep(NA, n),
b_uv = rep(NA, n),
b_yel = rep(NA, n),
sat_uv = rep(NA, n),
sat_yel = rep(NA, n),
ccar = rep(NA, n))
for(k in 1:nrow(color_met)){
color_met$b_tot[k] <- mean(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 701 &
pop_r2$wavelength > 319))
color_met$b_uv[k] <- mean(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 401 &
pop_r2$wavelength > 319))
color_met$b_yel[k] <- mean(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 626 &
pop_r2$wavelength > 549))
color_met$sat_uv[k] <- sum(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 401 & pop_r2$wavelength > 319)) /
sum(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 701 & pop_r2$wavelength > 319))
color_met$sat_yel[k] <- sum(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 626 & pop_r2$wavelength > 549)) /
sum(subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength < 701 & pop_r2$wavelength > 319))
color_met$ccar[k] <- (subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength == 700) -
subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength == 450)) /
subset(pop_r2$reflectance, pop_r2$id == color_met$id[k] &
pop_r2$wavelength == 700)
}
## Add little function for smoothed lines
my_fn <- function(data, mapping, method="loess", ...){
p <- ggplot(data = data, mapping = mapping) +
geom_point() +
geom_smooth(color = "coral3", fill = "coral3", alpha = 0.3, method=method, ...)
p
}
## Make a pairwise plot of relationships
p_cors <- GGally::ggpairs(color_met, columns = 2:7, lower = list(continuous = wrap(my_fn, method = "lm"))) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Make correlation heatmap
cormat <- round(cor(color_met[, 2:7]), 2)
cormat[lower.tri(cormat)] <- NA
melted_cormat <- reshape2::melt(cormat, na.rm = TRUE)
for(i in 1:nrow(melted_cormat)){
if(melted_cormat$Var1[i] == melted_cormat$Var2[i]){melted_cormat$value[i] <- NA}
}
melted_cormat <- subset(melted_cormat, is.na(melted_cormat$value) == FALSE)
p_cors2 <- ggplot(data = melted_cormat, aes(Var2, Var1, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "orange", high = "slateblue", mid = "white",
midpoint = 0, limit = c(-1, 1), space = "Lab",
name = "Person\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1),
axis.text.y = element_text(size = 12), legend.position = "top") +
coord_fixed() +
xlab("") + ylab("")
## Calculate variance along the x axis
var_pop <- data.frame(var = rep(NA, nrow(pop_r)))
for(i in 1:nrow(var_pop)){
temp <- t(pop_r[i, 1:(ncol(pop_r) - 1)])
var_pop$var[i] <- var(temp)
}
var_pop$wl <- lo_pred$wavelength
var_plot <- ggplot(var_pop, mapping = aes(x = wl, y = var)) +
geom_line(color = "coral3", lwd = 2) +
theme_bw() +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 12),
axis.title = element_text(size = 12)) +
xlab("Wavelength") +
ylab("Population variance")
## Put the output together
sim_out <- list(
wide_specs <- pop_r,
long_specs <- pop_r2,
color_mets <- color_met,
spec_plot <- p_spec,
cor_plot <- p_cors,
cor_plot2 <- p_cors2,
var_plot <- var_plot,
var_df <- var_pop
)
return(sim_out)
}Using the simulation
This is the easy part now that the function is written. We can compare the covariation patterns for a variety of different scenarios. Here I’m just changing the parameters in the functions and then showing the default plots that are produced by the function with everything else held the same. Note that because of the way the sources of variation are calculated sequentially and in some cases by multiplication, changing one parameter can still have some impact on other measures.
We could consider different mechanisms. For example, I’m sort of assuming here that a fixed amount of lutein works by absorbing a certain percentage of reflectance at each wavelength. An alternative would be that a fixed amount of luten absorbs a certain absolute amount of reflectance at each wavelength. I don’t really know which is a better assumption, but in some cases that will lead to different patterns. Similarly, in some cases it might matter which order these sources are added in, but I’ve always added them in the same order 1-4 here. For example, adding total reflectance changes after lutein might result in a different pattern. I don’t think that these choices will have a huge effect on the general conclusions, but something to think about.
The simulation is also a drastic oversimplification of the physical properties of the feather structures (i.e., they are not included at all). Because pigments have to be integrated with the feather matrix, it is totally possible that adding more pigment will change other characteristics of the shape (e.g., reducing total brightness) even outside of the range of absorbance of the pure pigment.
Finally, for the examples below I’ve just stuck with a single amount of variation added when an effect is ‘on’ (i.e., zero variation vs. about a 15% SD). Of course we could also do things like a small amount of variation in one source and large amount in another etc, at some point those would converge toward the single variation examples with the source providing the largest amount of variation having the biggest effect.
All of these examples are using the ‘after’ curve, but that choice doesn’t seem to make much difference.
1. Shape 1: boosting peak reflectance
shape1_sim <- carotenoid_sim(shape1_sd = 0.15)
grid.arrange(grobs = list(shape1_sim[[4]], shape1_sim[[6]], shape1_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
2. Shape 2: boosting UV reflectance
shape2_sim <- carotenoid_sim(shape2_sd = 0.15)
grid.arrange(grobs = list(shape2_sim[[4]], shape2_sim[[6]], shape2_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
3. Total reflectance with no shape change
A small difference here from the previous version is that I changed the order of the operations so that all of the white ‘base’ layer additions are covered first before the lutein subtraction is made. Previously I was shifting reflectance up and down after adding lutein. This approach is probably (?) a better approximation although we could discuss. The previous approach of just raising or lowering the entire curve without changing shape at all might be a better way to model measurement error…
totalb_sim <- carotenoid_sim(totalb_sd = 0.09)
grid.arrange(grobs = list(totalb_sim[[4]], totalb_sim[[6]], totalb_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
4. Carotenoid deposition
car_sim <- carotenoid_sim(pigment_sd = 0.15)
grid.arrange(grobs = list(car_sim[[4]], car_sim[[6]], car_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
Combinations
In the previous version I had a series of combinations of variation in 1/2/3/4 mechanisms all together. I can easily add those back in, but I think they were hard to gather much info from and I’ve taken them out for now.
Correlations
We discussed that one way we could recover some different patterns is if carotenoid deposition is correlated with feather structure/quality. The revised simulation can handle this. Any combination of covariance between the 3 or 4 mechanisms can be simulated, but I’m just starting simple here.
5. Total brightness covaries with carotenoids
Here I’m including a strong covariance between total underlying feather brightness and the amount of carotenoids deposited.
tb_car_sim <- carotenoid_sim(totalb_sd = 0.09, pigment_sd = 0.15, totalb_pigment_cor = 0.7)
grid.arrange(grobs = list(tb_car_sim[[4]], tb_car_sim[[6]], tb_car_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
6. Measurement error
No variation at all in any of above mechanisms, only measurement error that raises or lowers the entire spectra without changing shape.
error_sim <- carotenoid_sim(measure_error = 0.05)
grid.arrange(grobs = list(error_sim[[4]], error_sim[[6]], error_sim[[7]]),
widths = c(1, 1),
layout_matrix = rbind(c(1, 1),
c(2, 3)))
Comparing variance patterns
I’m pulling out the lines for population variance from the five simulations above for comparison. This is just a single comparison of 50 simulated individuals in each population. I could of course do more simulations or run in different ways, but thought it would be useful to see them together and to compare the correlation heatmaps side by side.
var_comp <- data.frame(var = c(shape1_sim[[8]]$var, shape2_sim[[8]]$var, car_sim[[8]]$var,
totalb_sim[[8]]$var, tb_car_sim[[8]]$var, error_sim[[8]]$var),
comp = rep(c("Shape1", "Shape2", "Carotenoid", "TotalB", "Tb_Car_Cor", "Error"), each = nrow(shape1_sim[[8]])),
wl = rep(shape1_sim[[8]]$wl, 6))
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
ggplot(var_comp, mapping = aes(x = wl, y = var, color = comp)) +
geom_line(lwd = 1.2) +
theme_bw() +
scale_color_manual(values = cbPalette) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title = element_text(size = 12)) +
ylab("Population variance") + xlab("Wavelength")
Comparison of correlation heatmaps from each simulation above
ggpubr::ggarrange(shape1_sim[[6]] + ggtitle("Shape1"),
shape2_sim[[6]] + ggtitle("Shape2"),
car_sim[[6]] + ggtitle("Carotenoid"),
totalb_sim[[6]] + ggtitle("Total B"),
tb_car_sim[[6]] + ggtitle("Tb ~ Carot"),
error_sim[[6]] + ggtitle("Error"),
nrow = 2, ncol = 3,
common.legend = TRUE)
Table comparing correlations
removed this for now since it is sort of replaced above
# I'm just compiling the rough output of all the combinations in the previous section into a table for comparison.
#
# - A = Variation in structural shape
# - B = Variation in UV shift
# - C = Variation in Lutein deposition
# - D = Variation in overall reflectance
#
# Note that I just added these by hand because I was having a hard timing comparing all the differences flipping back and forth so there could well be errors here. I don't have a strict guideline for what I called ++/w+/0/w-/--; I was just eyeballing the plots. It's also a randomization so the results will be slightly different each time. I've only run with 25 individuals.
# comps <- as.data.frame(matrix(nrow = 15, ncol = 14))
# colnames(comps) <- c("A", "B", "C", "D", "A_B", "A_C", "A_D", "B_C", "B_D", "C_D", "A_B_C", "A_B_D", "B_C_D", "A_B_C_D")
# rownames(comps) <- c("bt_v_buv", "bt_v_by", "bt_v_suv", "bt_v_sy", "bt_v_ccar",
# "buv_v_by", "buv_v_suv", "buv_v_sy", "buv_v_ccar",
# "by_v_suv", "by_v_sy", "by_v_ccar",
# "suv_v_sy", "suv_v_ccar",
# "sy_v_ccar")
# comps[, "A"] <- c("++", "++", "--", "++", "++", "++", "--", "++", "++", "--", "++", "++", "--", "--", "++")
# comps[, "B"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "++", "--", "--", "--", "--", "++")
# comps[, "C"] <- c("++", "++", "--", "--", "--", "++", "--", "--", "--", "--", "--", "--", "++", "++", "++")
# comps[, "D"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "++", "--", "--", "--", "--", "++")
# comps[, "A_B"] <- c("++", "++", "w+", "w-", "w-", "++", "++", "--", "--", "w+", "w-", "w-", "--", "--", "++")
# comps[, "A_C"] <- c("++", "++", "0", "0", "0", "++", "0", "0", "0", "0", "0", "0", "++", "++", "++")
# comps[, "A_D"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "++", "--", "--", "--", "--", "++")
# comps[, "B_C"] <- c("++", "++", "++", "--", "w-", "++", "++", "--", "0", "++", "--", "0", "w-", "0", "w+")
# comps[, "B_D"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "++", "--", "--", "--", "w-", "w+")
# comps[, "C_D"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "++", "--", "--", "--", "w-", "++")
# comps[, "A_B_C"] <- c("++", "++", "w+", "w-", "w-", "w+", "++", "--", "w-", "0", "0", "0", "--", "0", "++")
# comps[, "B_C_D"] <- c("++", "++", "++", "--", "--", "++", "++", "--", "--", "w+", "--", "--", "--", "--", "++")
# comps[, "A_B_D"] <- c("++", "++", "w+", "--", "--", "++", "++", "--", "--", "w+", "--", "--", "--", "w-", "++")
# comps[, "A_B_C_D"] <- c("++", "++", "w+", "--", "--", "++", "++", "--", "w-", "w+", "w-", "w-", "--", "w-", "++")
#
# comps %>%
# kbl() %>%
# kable_styling()